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Abstract 

We investigate the possibility of programming arbitrarily complex space-time patterns, and transitions 
between such patterns, by gene networks. We consider networks with two types of nodes. The w-nodes, 
called centers, are hyperconnected and interact one to another via M-nodes, called satellites. This central- 
ized architecture realizes a bow-tie scheme and possesses interesting properties. Namely, this organization 
creates feedback loops that are capable to generate any prescribed patterning dynamics, chaotic or pe- 
riodic, or stabilize a number of prescribed equilibrium states. We show that activation or silencing of a 
node can sharply switch the network dynamics, even if the activated or silenced node is weakly connected. 
Centralized networks can keep their flexibility, and still be protected against environmental noises. Find- 
ing an optimized network that is both robust and flexible is a computationally hard problem in general, 
but it becomes feasible when the number of satellites is large. In theoretical biology, this class of models 
can be used to implement the Driesch-Wolpert program, allowing to go from morphogcn gradients to 
multicellular organisms. 

1 Introduction 

The richness of Alan Turing's ideas hides somehow their unity. Is there a relation between the "chemical 
theory of morphogenesis" (Turing 1952) and the "universal machine" , or other, less known works, such as 
"Intelligent machinery" (Turing 1969, Teuscher and Sanchez 2001) in which he anticipates random binary 
networks? As emphasized by M.H.A. Newman (Newman 1955), a common denominator of Turing's 
scientific work is the quest for a mechanical explanation of nature. However, an even deeper unifying 
idea concerns the computability of nature and reciprocally, how nature computes. Turing looked for 
a mechanical support for natural pattern computation and found an analog machine, working by the 
chemical morphogens. Had he had known about gene networks, he would have probably analyzed the 
computational capacity of these networks to make patterns and multicellular organisms. In this paper we 
discuss a particular class of gene networks, the centralized or bow tie networks, and show that they can 
"compute" multicellular organisms and comply with important desiderata of life such as flexibility and 
robustness. 

Flexibility and robustness are important properties of living systems in general, most particularly 
observed during development from egg or embryo to a large, fully organized organism. Flexibility means 
the capacity to change, when environmental conditions vary. Opposite to this, robustness is the capacity 
to support homeostasis in spite of external changes. Intriguingly, biological systems are in the same time 
robust and flexible. Organization of the body plan in development, should be robust under unavoidable 
fluctuations of maternal gradients, embryo size, and environment conditions (for instance, temperature). 



This is a viability condition. On the other hand, developmental systems should be flexible in order to 
produce a number of different patterns and complicated dynamics. 

Pattern formation processes are important for the body plan establishment via cell fate decisions 
in multicellular organisms. Mathematical modeling of these phenomena (Turing 1952, Meinhardt 1982, 
Murray 1993, Wolpert et al. 2002, Mjolness et al. 1991, Reinitz and Sharp 1995, Page et al. 2005) is 
based on systems of reaction-diffusion equations, sometimes with spatial inhomogencous reaction terms. 
Two different approaches, both considering that laws of physics and chemistry are sufficient to account 
for making of a multicellular organism, are fundamental in modeling of body plan organization. The first 
approach, pioneered by the seminal work of Turing (1952), uses diffusion-driven instability as a patterning 
mechanism. For Turing's mechanism, a spatial dependence of the reaction term is not necessary, and the 
translation symmetry breaking needed for body plan organization results from the Turing instability. The 
second approach is based on Driesch- Wolpert positional information (Wolpert et al. 2002, Wolpert 1970). 
The corresponding models can be also based on reaction-diffusion equations, but in this case the models 
have space dependent reaction terms and no translation symmetry, therefore Turing instability is not 
needed. The main example of this type of models is the gene circuit model (Mjolness et al 1991, Reinitz 
and Sharp 1995). In this case, spatial organization is triggered by pre-patterns of signaling molecules, 
generically called morphogen gradients. It is remarkable that germs of this second approach can be found 
in the conclusion of Turing's 1952 paper, where it is suggested that "most of a organism, most of the time, 
is developing from one pattern into another, rather than from homogeneity into a pattern" . 

The body plan, considered as well defined, stable sequence of transitions from one pattern to another, 
can be encoded in a system of gene-gene interactions or gene network. For instance, in the segmentation 
along the anterior-posterior axis of Drosophila (fruit-fly) embryos, the chemical support of the pre-pattern 
is the maternal bicoid gradient developed in eggs soon after fecundation. This gradient induces spatially 
localized expression of segmentation genes (hunchback, kriippel, giant, knirps, tailless, fushi tarazu, even 
skipped, runt, hairy, odd skipped, paired, sloppy paired, etc.) forming a gene network. This gene network 
employs several types of interactions to stabilize the segmentation pattern. Some of these interactions 
originate in trans, i.e., far, on the DNA sequence, from the gene, and are due to transcription factors (TF) 
and microRNAs (miRNA). Other interactions originate in cis, i.e., close, on the DNA sequence, to the 
gene. Indeed, the zygotic genome contains cis-regulatory elements (CREM) controlling expression of the 
segmentation genes. The gene circuit model (Mjolness et al 1991, Reinitz and Sharp 1995) accounts for 
part of the trans interactions, considering that segmentation genes are mutually regulated transcription 
factors. The miRNA and the CREMs interactions are not represented in this model. Although the role in 
stabilizing the development has been experimentally proven for miRNAs(Li et al. 2009) and for CREMs 
(Ludwig et al. 2011), the mechanistic details of these interactions are still unknown. MiRNAs and CREMs 
can be abstractly considered as intermediate nodes in a gene network, mediating interactions between 
transcription factors. In such networks, TFs can be target hubs, being controlled by many CREMs, and 
also source hubs, because they can bind to many CREMs. Similarly, a few bioinformatics studies (Shalgi 
et al. 2007) suggest the existence of many genes submitted to extensive miRNA regulation with many 
TF among these target hubs. Direct testing of these interactions have recently shown that important 
transcription factors can be regulated by multiple miRNA (Tu and Bassler 2006, Martinez et al 2008, 
Wu et al 2010, Peter 2010). Without excluding other applications of our mathematical framework, we 
consider the TF-miRNAs networks as well as the TF-CREMs networks as possible examples of centralized 
networks, or bow-tie networks. 

The main goal of this paper is to show, by rigorous mathematical methods, the following new results: 

i centralized networks can create "a multicellular organism" consisting of many specialized cells where 
the network dynamics within each cell can have a different attractor, 

ii this pattern is robust under variations of morphogenetic fields; our system performs trade-offs be- 
tween flexibility and robustness, 

iii bifurcations between attractors can be obtained by gene silencing or reactivation. 

These results, however, would be useless, without algorithms that can resolve, in polynomial time 
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Poly{N), (where N is the gene number), the problem of a prescribed comphcated and robust pattern 
construction ( "computation of a robust organism" , CRO problem) . It is one of the key questions in 
development, why evolution had a sufficient time to construct complicated organs and organisms (Darwin, 
Origins of Species, Chapter 6). This problem is, in fact, a hard combinatorial one. Using new ideas in 
such problems, we show, under some assumptions, that 

iv for centralized networks with large hub connectivity, the CRO problem is feasible in polynomial 
Poly{N) time. 

Similar ideas, that bow-tie connectivity can play a role in flexibility and robustness, have been proposed 
by (Csete and Doyle 2004, Ma et al. 2007) in the context of metabolism, but lacked mathematical proofs. 
In theoretical computer science, it was shown that artificial neural networks can simulate any Turing 
machine (Sicgelmann and Sontag 1991, 1995). Also, it was shown that networks can simulate any time 
trajectories (Funahashi and Nakamura 1993) and any attractors (Vakulcnko 1994, 2000, Vakulcnko and 
Gordon 1998). 

We extend these results to simulation of any spatio-temporal structure, with any attractors. Since 
the pioneering ideas of Delbriick (1949), it became well accepted that differentiation and specialization 
of initially undifferentiated clone cells can be understood via multiple dynamical structures and attrac- 
tors (Thomas 1998). In particular, differences between gene expression programs can be understood as 
differences between attractors of dynamical gene networks. At least mathematically, the possibility to 
control any spatio-temporal pattern is equivalent to the possibility to organize any multicellular organism. 
Thus, we show that centralized networks can be used to implement Driesch-Wolpert positional information 
paradigm in order to organize a multicellular organism. This organism consists of a number of specialized 
cells, each cell type being dynamically characterized by distinct attractors. The complexity of the attrac- 
tors, that can be arbitrarily large, can be programmed by gradients of morphogens. Transitions between 
attractors can be performed by acting on key nodes of the network. Contrary to previous theories of 
random networks (Kauffman 1969, Aldana 2003, Aldana and Cluzel 2003), these key nodes do not have to 
be hubs. Furthermore, we show that patterning in such networks is maximally flexible in the sense that it 
can produce any structurally stable attractor. We also prove that (and show how) optimally flexible and 
robust structures can be computed in polynomial time and can thus be easily attained by evolution. 

The paper is organized as follows. Centralized networks are introduced in Section 2. A first theorem 
(Proposition 2.3) concerns with the flexibility of general centralized networks. We show that these networks 
are capable to generate practically all structurally stable prescribed dynamics, chaotic or periodic, and can 
have any number of equilibrium states. Another key result (Theorem 2.5) can be interpreted, in biological 
terms, as follows. The centralized networks are capable to create a "multicellular organism" , where each 
cell have a prescribed time dynamics. This assertion can be considered as a mathematical realization 
of the Wolpert approach since this intrinsic dynamics in a cell is predetermined only by the morphogen 
concentration in this cell. In Section 3 we show that gene activation or silencing can produce a sharp 
change of dynamics even if this gene is weakly connected in the network (it is well known that a mutation 
in a hub can sharply change the dynamics, see Aldana 2003). We show that in such a way one can obtain 
arbitrary bifurcations. In Section 4 we consider the robustness of centralized networks and show how these 
can acquire protection against environmental perturbations. We show that the design of a network that is 
both flexible and robust can be stated as an optimization problem for a discrete spin hamiltonian. When 
the number of satellites N is large, the optimization problem can be solved in polynomial time, Poly{N). 

2 Centralized networks 

By centralized networks we mean networks that contain a few strongly connected nodes (hubs) and a 
number of less connected, satellite nodes. A typical example is given by scale-free networks (Albert and 
Barabasi 2002, Lesne 2006), that occur in many areas, in economics, biology and sociology. In the scale- 
free networks the probability P{k) that a node is connected with k neighbors, has the asymptotics Ck~'^, 
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with 7 G (2, 3). Such networks typically contain a few hubs and a large number of satellite nodes. Hence, 
scale-free networks are, in a sense, centralized. 

In order to model dynamics of centralized networks we adapt a gene circuit model proposed to describe 
early stages of Drosophila (fruit-fly) morphogenesis (Mjolness et al. 1991, Reinitz and Sharp 1995). To 
take into account the two types of the nodes, we use distinct variables Vj, Ui for the centers and the 
satellites. The real matrix entry Aij defines the intensity of the action of a center node j on a satellite 
node i. This action can be either a repression Aij < or an activation Aij > 0. Similarly, the matrices 
B and C define the action of the centers on the satellites and the satellites on the centers, respectively. 
Let us assume that a satellite does not act directly on another satellite. We also assume that satellites 
respond more rapidly to perturbations and are more diffusive/mobile than the centers. 

Let M, N be positive integers, and let A, B and C be matrices of the sizes N x M, M x M and M x N 
respectively. We denote by Ai,Bj and Cj the rows of these matrices. To simplify formulas, we use the 
notation 

M A'l N 

AijVj = AiV, ^ Bjivi = BjV, ^ CjkUk = CjU. 

j = l 1=1 k=l 



Then, the gene circuit model reads: 

— - = diAi 

at 



+ fjcr [AiV + bim{x) - - X^Ui, (1) 



— djAvj + rjCj (BjU -f- CjU + bjinix) — hj) — XjVj, (2) 



dt 

where m{x) represents the maternal morphogen gradient, i — 1, N, j = 1, M . We assume that the 
diffusion coefficient di,di and maximal production rates ri^fi are non- negative: di,di,ri,fi > 0. Here 
the morphogenetic field to(x) and unknown gene concentrations Ui{x,t),Vj{x,t) are defined in a compact 
domain x € [dim{f2) < 3) having smooth boundary df2, x £ f2 and cr is a monotone and smooth (at 
least twice differentiable) "sigmoidal" function such that 



ct(-oo)=0, (t(+oo) = 1. (3) 



Typical examples can be given by 



^W- ,^ \ ,y a{h)^U^J= + l). (4) 
l + exp[-h) 2VV1 + /J / 

The function a{/3x) becomes a step-like function as its sharpness /? tends to oo. 
We also set the Neumann boundary conditions 

\7u,{x, t) ■ n{x) = 0, Vwj (cc, t) ■ n{x) = 0, {x e dH). (5) 

They mean that the flux of each reagent through the boundary is zero (here n denotes the unit normal 
vector towards the boundary df2 at the point x). Moreover, we set the initial conditions 

u,(a;,0) = ^^(a;) > 0, Vj(x, 0) = 0j(x) > {x e H). (6) 

It is natural to assume that all concentrations are non-negative at the initial point, and it is easy to show 
that they stay non-negative for all times (see below). 

Neglecting diffusion effects we obtain from (HJ,© the following shorted system: 

^ = ho- (^AiV + him{x) - - XiUi, (7) 
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= r-jcr [BjV + Cju + bjm{x) - hj) - XjVj. (8) 

This is a Hopfield-likc network model (Hopfield 1982) with thresholds depending on x (contrary to the 
Hopfield model, the interaction matrices are not necessarily symmetric). In this case we remove all 
boundary conditions ([5]). If only di = we remove the corresponding boundary conditions for Vi. 



2.1 Existence of solutions 

Let us introduce some special functional spaces (Henry, 1981). Let us denote H = L2{fi)'"' the Hilbert 
space of the vector value functions w. This space is enabled by the standard L2- norm defined by 
\\w\\^ = J^\w{x)\^dx, where \w\^ = X)^!- ^or a > we denote by Ha the space consisting of all 
functions w £ H such that the norm ||w||q is bounded, here ||w||^ = ||(— ^ + /)"w|p. These spaces have 
been well studied (see Henry 1981 and references therein). The phase space of our system is H = {w = 
(u,v) : u £ V £ H}, the corresponding natural fractional spaces are denoted by Ha and Ha, here 
Hq = H and Hq = H. Denote by Ba{R) the n-dimensional ball in Ha centered at the origin with the 
radius R: Ba{R) = {w : w € Ha, \\w\\a < R}- 

In our case all fi(w,x) are smooth in w,x, therefore, the standard technique (Henry 1981) shows that 
solutions of ([3]) exist locally in time and are unique. In fact, our system can be rewritten as an 
evolution equation of the form 

wt = Aw + f{w), (9) 

where / is a uniformly bounded map from Ha to H (since sup^.^^ \w\ < c\\w\\a, a > 3/4, and the 
derivative cr'(z) is uniformly bounded in z) and a linear self-adjoint negatively defined operator A generates 
a semigroup satisfying the estimate || exp(Af)u'|| < exp(— /3t)||w|| with a /3 > 0. 

Let us prove that the gene network dynamics is correctly defined for all t and solutions are non-negative 
and bounded. In fact, there exists an absorbing set B defined by 

B = {w ^ {u,v) -.0 < Vj < rjXj\ 0<u,< fi~Xi\ j = 1, ...,M, i = 1, ...,N}. 

One can show, by super and subsolutions, that 

0<Ui{x,t) < 0i(a;)exp(-Ait) + fiA,"^(l - exp(-Ait)), ^^^^ 
< v^{x,t) < (j>i{x)ejip{-Xit) + nX^^l - exp(-A,t)). 

Therefore, solutions of (IT|), Q not only exist for all times t but also they enter the set B at a time moment 
to and then they stay in this set for all t > to. So, our system defines a global dissipative semiflow (Henry, 
1981). 



2.2 Reduced dynamics 

The key idea is to find a simpler asymptotic description of system dynamics. It is possible under some 
assumptions, we suppose here that the w-variables are fast and the u-ones are slow. We show then that the 
fast u variables are captured, for large times, by the slow v modes. More precisely, one has u = U{v) + it, 
where {t is a small correction. This means that, for large times, the satellite dynamics is defined almost 
completely by the center dynamics. 

To realize this approach, let us assume that the parameters of the system satisfy the following condi- 
tions: 

\A,il\B.al\a,l\U\h,\<Co, (11) 
where i = 1, 2, N, ij = 1, M, j = 1, N, 

< Ci < Aj, dj < C2, (12) 
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\bj\,\h\<C3, sup|m(x)| <C4, (13) 

and 

n = nRi, fi = nRi, (14) 

where 

<C5, \\\<Cq, (15) 

dj = kJj, <dj < Cy, (16) 

where k is a small parameter, and where all positive constants Ck are independent of k. 

Proposition 2.1. Assume the space dimension Dimfi < 3. Under assumptions ill]). ^14^ 
for sufficiently small k < kq solutions (it, v) of (0), (0), and (0) satisfy 

u = U{x,v{x,t)) + u{x,t), (17) 
where the j-th component Uj of U is defined as a unique solution of the equation 

djAUj - Xj Uj ^nGjiv), (18) 
under the boundary conditions where 

Gj = Rja (^Ajv{x, t) + bjm{x) — hj 

The function u satisfies the estimates 

||w|| + ||Vw|| < ciK^ +i?exp(-/3t), /3>0. (19) 
The V dynamics for large times t > Ci \ logKj takes the form 

dv ' 

--} = KF,iu,v) + w,, (20) 
ot 

where Wi satisfy 

\\Wi\\ < CqK^ 

and 

Fi{u,v) = diAvt + i?icr(B,w + CiU{x,v) + b^m - hi) - \,Vi. 

Constants co,Ci do not depend on k as k ^ but they may depend on Ri,Ri,Ci. 

A tedious proof of this assertion is basically straightforward; it is based on well known results (Henry 
1981) and is relegated to the Appendix. 

An analogous assertion holds for shorted system (l7]),([8]). In this case the functions Ui can be found by 
an explicit formula. Namely, one has 

U^{x,v{x,t)) = kV^, = R,~X:[^a (^Ajv{x,t) +bjm{x) -h^Y (21) 

For large times the reduced v dynamics has the same form ((20|) with d,; = 0. 
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2.3 Realization of prescribed dynamics by networks 

Our next goal is to show that dynamics (j20[) can realize, in a sense, arbitrary structurally stable dynamics 
of the centers. To precise this, let us describe the method of realization of the vector fields for dissipa- 
tive systems (proposed by Polacik (1991), for applications see, for example, (Dancer and Polacik 1999, 
Rybakowski 1994, Vakulenko 2000). One can show that some systems possess the following properties: 

A These systems generate global semifiows Sj, in an ambient Hilbert or Banach phase space H . These 
semiHows depend on some parameters V (which could be elements of another Banach space B). They 
have global attractors and Enite dimensional local attracting invariant - manifolds M. , at least for 
some V. 

(Remark: in some cases, these manifolds can be even globally attracting, i.e., inertial. Theory of 
invariant and inertial manifold is well developed, see (Marion 1989, Mane 1977, Constantin et al 1989, 
Chow and Lu 1988, Babin and Vishik 1988). 

B Dynamics of Sy, reduced on these invariant manifolds is, in a sense, "almost completely controllable" . 
It can be described as follows. Assume the differential equations 

f^=F{p), FeC\B") (22) 

define a dynamical system in the unit ball B" C R". 

For any prescribed dynamics iT^) and any 6 > 0, we can choose suitable parameters V = V{n,F,5) 
such that 

Bl The semiflow has a C^- smooth locally attracting invariant manifold Ai-p diffcomorphic to i?"; 
B2 The reduced dynamics S^y,\Mv defined by equations 

^=%,P), FeCi(B") (23) 

where the estimate 

\F - F\cHB^^) < S (24) 

holds. In other words, one can say that, by V , the inertial dynamics can be specified to within an 
arbitrarily small error. 

Thus, all robust dynamics (stable under small perturbations) can occur as inertial forms of these 
systems. Such systems can be named maximally dynamically flexible, or, for brevity, MDF systems. 

Such structurally stable dynamics can be "chaotic" . There is a rather wide variation in different 
definitions of "chaos" . In principle, one can use here any concept of chaos, provided that this is stable 
under small -perturbations. To fix ideas, we shall use here, following classical tradition (Ruelle and 
Takens 1971, Newhouse, Ruelle and Takens 1971, Smale 1980, Anosov 1995), such a definition. We say 
that a finite dimensional dynamics is chaotic if it generates a hyperbolic invariant set 7^, which is not 
a periodic cycle or a rest point. For a definition of hyperbolic sets see, for example, (Ruelle 1989); a 
famous example is given by the Smale horseshoe. If, moreover, this set F is attracting we say that F 
is a chaotic (strange) attractor. In this paper, we use only the following basic property of hyperbolic 
sets, so-called Persistence (Ruelle 1989, Anosov 1995). This means that the hyperbolic sets are, in a 
sense, stable(robust): if generates the hyperbolic set F and S is sufficiently small, then dynamics ([25)1 
also generates another hyperbolic set F. Dynamics (22) and restricted to 7^ and F respectively, are 
topologically orbitally equivalent (on definition of this equivalence, see Ruelle 1989, Anosov 1995). 

Thus, any kind of the chaotic hyperbolic sets can occur in the dynamics of the MDF systems, for 
example, the Smale horseshoes, Anosov flows, the Ruelle- Takens-Newhouse chaos, see (Newhouse, Ruelle 
and Takens 1971, Smale 1980, Ruelle 1989). Examples of systems satisfying these properties can be 
given by some reaction diffusion systems (Dancer and Polacik 1999, Rybakowski 1994, Vakulenko 2000). 
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Although not yet observed in gene networks, structurahy stable chaotic itineracy is thought to play a 
functional role in neuroscience (Rabinovitch 1998). 

Let us apply this approach to network dynamics using the results of the previous section. To this end, 
assume that [T5)) and (fTO)) hold. Moreover, let us assume 

bi = Kbi, hi = ^^h^ (25) 

\i = K^Ai, di = K^di (26) 

where all coefficients hi and hi arc uniform in k as k — 0. These assumptions are useful for technical 
reasons. We also assume that all direct interactions between centers are absent, B = 0. This constraint 
is not essential but facilitates notation and calculations. 

Since Uj = 0{k,) for small k, we can use the Taylor expansion for cr in (|20p . Then these equations 
reduce to 

— '—^ — = diAvi + pi{CiV{x, v) + bim{x) - hi) - \iVi + Wi{x, t), (27) 

OT 

where Pi{x) = ficr'(O), i = 1, 2, M and r is a slow rescaling time: r = n^t. Due to conditions (|25p and 
(PS)) corrections Wi satisfy 

WwiW < CK. 

Let us focus now our attention to non-perturbed equation (j27p with Wi = 0. Let us fix the number of 
centers M . The number of satellites N will be considered as a parameter. 

The next important lemma follows from known approximation theorems of multilayered network theory, 
see, for example, (Barron 1993, Funaliashi and Nakamura 1993). 

Lemma 2.2. Given a number 5 > Q, an integer M and a vector field F = (i^i, ...,Fm) defined on 
the ball B^^ = < 1}, Fi e C^{B^^), there are a number N, a N x M matrix A, a M x N matrix C 
and coefficients hi, where i — l,2,...,Af, such that 

\F,{-)~C,W{-)\c^iBM)<5, (28) 

where 

Wi{v) ^ a {A,v - hi) , (29) 

where v = (wi, vm) G R*^. 

This lemma gives us a tool to control network dynamics and patterns. First we consider the case when 
the morphogens are absent. Formally, we can set bi = bj = di = 0. Assume hi = 0. Then equations ^7} 
with Wi ~ reduce to the Hopfield-like equations for variables Vi = Uj;(t) that depend only on r: 

^=KiWiv)-~Xm, (30) 
dr 

where I — 1, A/, the matrix K is defined by Kij ~ piCijRjXj^ . The parameters V of ([50)1 are K, Af, 
hj and Xj. 

In this case one can formulate the following result. 

Proposition 2.3. Let us consider a -smooth vector field Q{p) defined on a ball _B^^ C R*^ and 
directed strictly inside this ball at the boundary dB^ : 

F{p)-p<0, pedB^^. (31) 

Then, for each S > 0, there is a choice of parameters V such that fi30(l 5 -realizes the system 112!^) . This 
means that is a MDF system. 

This proposition follows from the Prop. 2.1 and Lemma 2.2. 

Prop. 2.3 implies the following important corollary: all structurally stable dynamics, including periodic 
and chaotic dynamics can be realized by centralized networks. The proof of this fact uses the classical 
results on the persistence of hyperbolic sets, and on the existence of invariant manifolds (Ruelle 1989), 
see (Vakulenko 2000). 
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2.4 Pattern and attractor control by Wolpert positional information 

Above we have considered a spatially homogeneous case. Proposition 2.3 shows that a centralized network 
can approximate an arbitrary prescribed dynamics. Thus, it is shown that cells can be programmed to 
have arbitrarily complex dynamics. By network rewiring or by interaction tuning, one can switch between 
various types of dynamics. During development these switches are position dependent, and induce cell 
differentiation into specific spatial arrangements. 

Let us show that the centralized networks, coupled to morphogen gradients, can generate any spatio- 
temporal pattern as support for multicellular organization. We consider shorted dynamics ([7]), (O that 
is reasonable for cellularized developmental stages, where cell walls prevent a free diffusion of regulatory 
molecules. Although other phenomena such as cell signalling can also lead to cell coupling, we do not 
discuss these effects here. 

Assume cell positions are centered at the points x € X = {a;i,X2, ...^Xk}, dimQ = 1, A" is a discrete 
subset of [0, L\. Let us show that eqs. dZ])-® can realize different dynamics at different points xi of the 
domain Q = [0, L]. 

We can formulate now the following. 

Theorem 2.5. (On translation of positional information into complex and variegated cell 
dynamics, or programming of multicellular organism). Suppose x E [0,L] C R and m{x) is a 
strictly monotone smooth function. 

Assume that < xi < X2 < ■■■ < Xk < L and that F^''^{p), I = l,2,...,fc is a family of -smooth 
vector fields defined on a unit hall C . We assume that each field defines a dynamical system, 
I.e., FW are directed inwards on the boundary dB . 

Then, for each S > there is a parameter V choice such that for shorted dynamics one has 

u = U{x,v) + u, 

where 

\u\ < Cexp(-/3r) + cti^. 
For x = Xi and for sufficiently large times the dynamics for v{xi,t) can be reduced to the form 

^^mxi,p), (32) 

where 

sup \Fixi,p)-F'^''>ip)\<S. (33) 

Here pi (r) can be expressed in a linear way via Vi ixi , r) by 

Vi{xi,T) - hm{xi) ^ paPtir). 

This theorem can be considered as a mathematical formalization of positional information ideas. It 
extends Driesch- Wolpert theory by incorporating gene networks and coping with their information process- 
ing role. Flexible gene networks have different dynamics and attractors, for different local concentrations 
of morphogens. The attractor selection ensures the cell fate decision. Concerning the relation between 
attractors and cell fate determination, see (Delbriick 1949, Thomas 1998). 

To prove this assertion, let us turn to eqs. (|27|) . where, taking into account biological arguments given 
above, we set di — 0,bi — 0. Let us set, to simplify formulas, pj ~ 1, hj ~ and Xj ~ 1. Then 

V,{v)=RfX-'a{A,v~h,). 
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Denote by Qi the sums Qi{v) ~ X^jli ^ij^ji"^) — ^iV- Removing the terms Wi in (P7)). one obtams 
that eqs. ([27]) reduce to 

= Q,{v{x,T))+h,m{^)-v^{x,T). (34) 
Let us fix a X = xi € X . Let us make the substitution Vi{xi, r) = ^^(t) + hm{xi) in p4p that gives 



i{z + bm{xi)) ~ (35) 



where & = (&i, &Af ), « = 1, M. 

Now we again use approximation Lemma 2.2. Let us consider a family of vector fields C^-smooth 
vector field F defined on a unit ball B^^ = {z € R^^, \z\ < 1} and directed strictly inside this ball at the 
boundary dB^^ : 

F(')(z)-z<0, zedB^^. (36) 

Assume m{x) is a strictly monotone function in x. The main idea is as follows: since all m[xi) = 
and rn{xj) = fij are different for j =^ I, the vector fields Q^^\z) = (5(z + bfxi) can approximate different 
vector fields _F''-'(z) for Z = 1, k and for z such that |z| < po, where po = 5 inini,j,/j-^; |6i||yu(xj) — fj.{xi)\. 

For each e > we can find an approximation Q satisfying 

sup \Q{z + bm{xi))~{poF^^'>Po^z) + z)\<pQe, (37) 

|z|<Po 

and 

sup \V{Q{z + bm{xi)))-^VpoF<-'\p^^z) + z)\< poe. (38) 

|2|<P0 

Then equation ((35)l reduces to 

J^PoFW(p„-1z) + Po.F«(z), 

where 

sup |F^')(z)| < 1, sup |VF(')(z)| < L 

^ =F«(p) + eF«(p). (39) 
ctr 

Let us notice that, if e is small enough, then for each index I, due to assumption p6[) . the trajectory 
p(i,p(0)) stays in when the starting point lies in B^^ : p{0) S i?^^. Consequently, our approximations 
(|57| give vector fields that, by ([M]). realize different dynamics for each a;;. 



We set Zi = poPi- This gives 



3 Sharp genetic switch by sateUite silencing/reactivation 

In the context of scale-free random networks, it was proposed (Aldana 2003) that removing of a strong 
connected center can sharply change the network attractor. Here we will show that one can obtain 
transitions between all possible structurally stable attractors by a single event acting on a specially chosen 
weakly connected satellite. Such a satellite interacts only to one or two centers. Such event may be, either 
deletion, silencing, or reactivation. Therefore, such a node can serve as a switch between two kinds of 
network behavior. Each of the type of behavior can be defined, for example, by an attractor or several 
coexisting attractors that can be fixed points, periodic or chaotic attractors. 
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To formalize these ideas mathematically, let us consider a system of ordinary differential equations 

^ = F{p,s), peB"cR" (40) 
at 

depending on a real parameter s. Here p = (pi,p2, ■■■,Pn), -B" is the unit ball centered at p = 0, and F is 
C^-smooth vector field directed inside the ball at the ball boundary for each s (sec ([M)) ). Let us consider 
sq,si such that sq ^ si and suppose that (j40[l has different attractors Aq and Ai for s = sq, s ~ si 
respectively. 

Consider, for simplicity, the gene circuit model ([1]), ^ without diffusion and space variables: 

^ = fjCT (^AiV - hi^ - XiUi, (41) 
dv ' 

=rj<j{CjU^hj)-X,Vj, (42) 

where i = 1, M+ 1, j = 1, N . The parameters V of this system are M, N, hi, hi, Xi,ri,fj,Xj and the 
matrices A, C. We can assume, without loss of generality, that we eliminate (by silencing) the M + 1-th 
satellite node, i = M + 1. As a result of this elimination, we obtain a similar system with i = 1, A/ 
and shorted matrices A, C. Of course, we can also consider the opposite event, which is to reactivate the 
M + 1-th node and recover the initial system this way. 

Theorem 3.1 For each e > there is a choice of the parameters V such that system {^1% with 
M satellite nodes e -realizes with s = sq and system j/^ j| ), with M + 1 satellite nodes e -realizes 
(|/^0| ) with s ~ si. 

To prove it, we use the following extended system 

^^=pF{p,s), peB" (43) 
ds 

f[s,(3)-i,s, SGR (44) 

where v > 1 and /(s) is a smooth function, /3, p > arc parameters. Then equilibrium points s^q of (|44p 
are solutions of 

f{s,P)^ys (45) 

The point Seq is a local attractor if fg{seq) < v- Let us denote Seq{Pk) = Sk, where fc = 0, 1, and let these 
roots of ps)) be stable, i.e., /'(sfe) < i^. 

Then, if p > is small enough, and Sk is a single stable rest point, the fast variable s approaches at 
SeqiP) and for large times t the dynamics of our system (^5]) . is defined by the reduced equations 

^=pFip,s,qm. (46) 

Now let us set 

/(s,/3) =2/3V(fe(s-M), ho<0 (47) 

where 6 is a large parameter. Then / is close to a step function with the step 2/3^. Therefore for Seg(/3) 
one has the asymptotics s^q = 2/3^i/~^ 4- 0(exp(— 6)) as 6 — oo. Thus, we can adjust parameters (3,b > 
in such a way that (|45p has a single stable root sq and the equation 

f{s,P/V2) = us (48) 

also has a single root si = /3'^ /ly ^ so- 
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Dynamics (|33]), (|33|) with / from (PT]) can be realized by a network (jH]), (|^^ in such a way. We 
decompose aU satelhtes Ui into two subsets. The first set contains satehites ui,U2, ...,um-i, the second 
one consists of the satehites um,um+i ( to single out this variables, let us denote um = yi,UM+i =2/2)- 
The main idea of this decomposition is as follows. We can linearize equations for the centers Vj assuming 
that the matrix C is small and B = (as above in Section 2). 

The y satellites realizes the dynamics (H^ by a center s: 

ds 

— ^^i,s + P{yi+y2), (49) 

^ = ~yk+Mb{s~ho)), fc-1,2. (50) 

Here we assume that i^, j3 is small enough, therefore, for large times this system reduces to (|44|l with / 
defined by (|47|) . We sec that this dynamics bifurcates into (|44)) with / = 0^a{h{s — Hq)) if wc remove 2/2 
in the right hand side of 

The rest of the equations, after a notation modification and linearization, take the following form 

^ ^ f,a (^A,v + D,s - h,^ - \,u„ i = l,...,A/-l (51) 
du ■ 

-^ = ~^3^j+rjCjU~h,, (52) 

where i = 1, M + 1, j = 1, N. 

Equations ((5T|) . ((52|) can e-rcalizc arbitrary systems (|40)) with the parameter s which can be shown as 
above (see Section 2), and this completes the proof. 



4 Robust dynamics 

Our definition of robust dynamics is inspired from similar ideas in viability theory (Aubin et al. 2005). Let 
us suppose that the dynamics (the global semiflow S**), generates a number of attractors. Each attractor 
A has an attraction basin B(A), that is an open set in the phase space. Assume that our initial data (f> 
lie in an attractor, (f> € A, and let us add some noise ^ to the dynamics, representing the effect of the 
environment. Trajectories become random, and then it is possible that, under this noise, the trajectory 
leaves B{A). 

We can now define the following characteristic of stability under the noise. Let us denote P{T, B{A), (f) 
the probability that the trajectory u(t, (j)) such that ii(0) = (j) ^ A stays in B{A) within the time interval 

[o,r]. 

Definition. Let us consider a network dynamics depending on some parameters V and on the noise 
^. We say that the dynamics is robust under the noise ^, if for each T and S > there is a choice of the 
parameters such that 

P{T,B{A),cp)<5 

for each attractor A and (j) ^ A. 



4.1 Centralized motif with noise 

For the rest of this section we consider a simplified network, with a single central node interacting with 
many satellites. This motif can appear as a subnetwork in a larger centralized network. In order to study 
robustness, we consider the case when the satellites and the center are under the influence of noise. More 
general situations, including perturbations of several centers and satellites, will be studied elsewhere. 
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The network dynamics is described by the following equations: 

= diAui - XiUi + a{biV - hi + ^i{x,t)), i^l,...,N, (53) 

dv ^ 

— = doAv - Xqv + cr(^ ttiU^ -ho + ^q{x, i)), (54) 

i—l 

The random fields £,i{x,t) summarize the effect of various extrinsic noise sources. These can be random 
variations of tlic morphogen, or environment noise, or genetic variability affecting network interactions. 

Intrinsic noise, resulting from stochastic gene expression, could be represented as supplementary terms 
outside the sigmoid function. In order to avoid further some tedious technical difBculties, we postpone 
the discussion of intrinsic noise to future work. Some aspects of the robustness of patterns with respect 
to intrinsic noise was studied numerically by Scott et al. (2010). 

The flexibility of such simple networks results from the preceding section. 

We can formulate the following problem: how to choose a network motif, robust under a given environ- 
mental noise, and simultaneously flexible? The choice can result either from genetic changes (for instance 
mutations, deletions or duplications of DNA regions) or from network plasticity (epigenetic changes, such 
as methylation and chromatin remodeling). 

Assume that the considered process is a choice of satellites Ui from a large pool of possible regulators. 
We can present this process as a choice of n indices = 1, n from a larger set In = {1,2,..., N} of 
indices, where N > n. This choice can be done by boolean variables Sj that multiply the coefficients aj: 
the j-th reagent participates in the network if sj = 1 and does not participate if s.; = 0. Let us make an 
important assumption allowing us to obtain a thermodynamical limit as N oo. We assume that 

\a,\<CN-\ N^oo. (55) 

Now we transform eqs. ([5^ . ([5^ . using the results of the previous subsection. Let v = q{x) and 
Ui = Ui{x) be equilibrium solutions of ([5 ^ . ([M)) where the noises are removed. We suppose that the 
assumptions of the previous subsection hold. Let us set 

V = q + v, = J7i + Ui, 

and U,u denote vectors [Ui, ...,Um), (wi, wat). Let us set temporarily Co = ( below we shall show 
how one can stabilize the system state, when 7^ 0). This gives 

^ = doAi - Xov + a{p{U + u) -h)~ a{p{U) - h), (56) 

where we use, for brevity, the notation p{u) — X^iLi SiUiUi. 
The second part of equations takes then the form 

= -diAui - XiUi + a{hi{q + u) + - hi) - a{hiq - hi). (57) 

To investigate equations ([55)1 . ([57|) . we use a special method justified in a rigorous way in Appendix. 
This holds under the following assumption: 

Assumption 4.4. The "morphogenetic" noises £,i{x,t) are independent on t: 

Ux,t)=i,{x), ^ = 0,1....,A^. 
The functions arc continuous in x and a priori bounded 

sup|e4(x)| <a, i = 0,l,...,iV. (58) 
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where a positive constant C* may be large but it is independent of A'^ for large N. 

Notice that Assumption 4.4 guarantees global existence of solutions Ui{x,t),v{x,t) of eqs. (j56[) . ((57)) 
for all t > 0. 

Intuitively, one can expect that the term v in (|57p in a is small and can be, thus, removed. Following 
this idea, let us introduce r]i as solutions of 

= -diArji - Ai?/, + a{hiq + C» - /^j) " (^{hq - hi). (59) 

If S,i are independent of t, and sufficiently regular in x then arguments of the previous section show that 
for large t 

n,{x,t)^f],{x), (60) 

where fji are solutions of elliptic equations 

d,Afi^ + Kn, = G,(e,(a;)), G.(e.(a:)) = cr(fo»g + - /».) - - (61) 

under zero Neumann boundary conditions. 

Let us consider equation (|56p for w. Assume p is small. Then we can linearize the nonlinear contribu- 
tions in the right hand side of this equation: 

a{p{U + u)-h)- (j{p{U) -h) = (t'{U - h)p{u) + 0{p{uY). 

We assume that Ui are close to fji. Thus, p{u) ~ p(j])- Calculations presented in the Appendix show that 
the fluctuation influence can be estimated through the quantity 

(5(s,T)= sup H(s,0, U{s,t) = \\p{f])\\^. (62) 

tG[0,T] 

If are independent of t, for large t one has 

His,t)^Uis) = \\pim'- (63) 
Notice that H can be rewritten in the form 

N N 

= II II W^,n'{fj{-))s^s,^,, (64) 

m— 1 m' — l 

where M^mm are random and 

{Vm, Vm'), (65) 

here (/, g) denotes the inner scalar product in H: (/, g) = J^^ fgdx, where dx is the standard Lebesgue 
measure. 



4.2 Hard combinatorial problems in network evolution 

We assume that Assumption 4.4 holds. The minimization of H(s) with respect to s should be done under 
the condition that at least one satellite is involved, i.e., 

N 

i?o(s) =iV-i^s, > 0. (66) 

i=l 

The analysis of the minimization problem for this random Hamiltonian is a computationally hard problem 
advanced firstly by methods from statistical physics of spin glasses (see, for example, (Mezard Zecchina 
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2002) for applications to hard combinatorial problems, and (Talagrand 2003) for rigorous justification). 
To make the analogy with spin glasses more transparent, we can make change Si = 2Si + 1, where spin 
variables Si take values 1 or —1. However, our problem is even more complicated because, besides (|66p. 
some other restrictions should be taken into account. 

In addition to (|66p , we must take into account restrictions connected with generation of several steady 
states qi, q2, Qm, to provide flexibility. Let us take a small e > 0. By adjusting Si we would like to 
obtain a set of equilibria close to qi . This gives the following restrictions 

N 

sup cr{y^ Sia^<7{biqi ~ h) ~ Xoqi - ho) < e, I ^1,2,...,M (67) 
.en ^ 

or, in a simpler form, 

sup \Ri{s, x) - Bi{x)\ < ce, / = 1,2,...,M (68) 

where 

N 

Ri = y^^MiiSj, 

Mii = aia{biqi - h), Bi = cr"'^ {X^qi - ho). 

Although Ri are linear in Si functions, the left hand side of (|68p is, in general, a nonlinear function of 
a complicated form. To overcome this difficulty, we replace the sup in (j68p by the L2- norm that gives 
quadratic in s functionals: 

Riis)^\\Riis,x)^Bi{x)\\^ <ce\ l^l,2,...,M (69) 

We use Lagrange multipliers /3; to take into account conditions ((69|) . This leads to the following 
Lagrange function: 

L 

F(s) -H(s)+^AR,(s)^ (70) 

/=i 

Let us remind that the matrix Wmm', that determines our hamiltonian H, is a random matrix depending 
on random fields ^i(x) through p{f]). Let us consider these fields as elements of the Banach space C°(i7) 
of all bounded continuous in x vector valued functions. Let ji^ be a probability measure defined on the 
subset of all such functions satisfying ([55)) . 

We also propose that variables s are chosen by a stochastic algorithm. The stochastic algorithm 
depends on some set of parameters P that can be adjusted. Let be a probability measure associated 
with this algorithm (this measure is defined below). 

We would like to have a small value of F for a "most part" of field ^ and s values, with respect to the 
product measure ^ = fi^ x fxp. 

Finally, the combinatorial problem can be formulated as follows: for a small number 5, find parameters 
P such that the probability (computed by the measure fi), 

Prob{Fia-),s)>S} (71) 

is small enough. 
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4.3 Mean field solution can be obtained by quadratic optimization 

We show here that the optimization problem is feasible when N is large. To this end, we define the mean 
field Lagrange function F that is obtained from F(^(-),s) by averaging with respect to fi. In order to 
estimate the deviations of F from F we use the Chcbyshev inequality: 

P{F, s) = Prob{\F{^{-),s) - F(.s)| > a} < a^'^VarF, (72) 

where the probability, the average and the variance should be computed by fi. 

The stochastic algorithm for choosing the satellites can be a simple Bernoulli scheme. Namely, let us 
consider Si as mutually independent random variables such that 

Pro6{sj = 1} = Pi. 

Thus, the mean field Lagrange function reads 

N N L 
i=l j=l 1=1 

where Wij is obtained from Wij by averaging with respect to fi^. Our main idea is as follows. 
Step 1: Quadratic programming for the mean field Lagrange function 

First, we minimize F with respect to pi. This is a quadratic programming problem that can be solved 
in polynomial time. 

QP to find a minimum F(p) under conditions 

< K < 1, (74) 

N 

^o(p)=^P, >0. (75) 

4=1 

The last condition is trivial and can be omitted. Therefore, we look for a minimum of a positively 
defined quadratic form on the multidimensional box. The well-known L.Khachiyan ellipsoid algorithm for 
this problem runs in Poly{N) time. This proves such a lemma: 

Lemma 4.5. // a solution of the problem QP exists, then it can he found in Poly{N) time. 

Step 2: Obtain a small variance of the Lagrange function in the limit N large 
Let us suppose that F < 5/2. Then, using ([72]) we get 

Pro&{F(^(-), s) >5] < Prob{\F{i{-),s) - F(s)| > 5/2} < 4:5-'^VarF. (76) 

Now let us estimate VarF. We consider VarH, the rest terms R; can be considered in a similar way. 
First we estimate variation with respect to s by the measure fip. One notices that 

VarH = E ^ SiSjSi>Sj>WijWi'j> — E'y^^SjSjWjjE^y^^ Sj'Sj'Wj'j' . 

iji'j' ij i'j' 

Notice that ii i i' and j ^ j' then 

Moreover, \ Wij \ = 0{N~'^) due to our assumption (|55p on and Assumption 4.4. Thus we have maximum 
of non-zero terms in DH , which have the order 0{N~^). Thus, the complete variation satisfies 

VarFl < CoN-\ (77) 
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where Co is uniform in as — > oo. 

Thus, for large N one has VarF — > 0, thus the probabihty (|7ip is arbitrarily small. This shows that 
the problem of minimization ([7^ is feasible in polynomial time Poly{N), when N is large enough. More 
precisely, we have the following 

Proposition 4.6. // a solution of the problem QP exists and N is large enough, then a solution s 
satisfying all restrictions and minimizing F at level 6 with a probability, arbitrarily close to I, can be found 
in Poly{N) time. 

Remark. Above we have studied the case = 0. For smooth S,q{x) we can obtain a robustness with 
respect to variations in a simple way. Namely, for large A^ one can choose the constant C in ([55]) ) large 
enough, then X^ili ~ >> l?o(a;)|. 

There arises, however, a natural question: how genetic networks can realize these sophisticated algo- 
rithms which arc capable to optimize the network robustness? A possible answer to this question is that 
Si could be themselves involved in a gene network of the form ([T]), We showed that gene networks 
are capable to simulate all structurally stable dynamics. The fact that this is equivalent to simulating 
arbitrary Turing machines and thus arbitrary algorithms follows from results of Koiran and Moore (1999). 

5 Conclusion 

We are concerned with dynamical properties of networks with two types of nodes. The w-nodes, called 
centers, are hyperconnected and interact one to another via many u-nodes, called satellites. We show, by 
rigorous mathematical methods, that this centralized architecture, widespread in gene networks, allow to 
realize two fundamental biological strategies: flexible and robust bow-tie control and Wolpert positional 
information concepts. 

We show how a combination of these strategies leads to the remarkable possibility to create a "mul- 
ticellular organism", where each "cell" can exhibit a complicated time behaviour, different for different 
cells. Centralized network architectures provide the flexibility important in developmental processes and 
for adaptive functions. 

Contrary to previous works on centralized boolean networks (Aldana 2003), we show that arbitrary 
bifurcations between attractors can be controlled by action on satellites, instead of actions on centers. 

To check the robustness of such architectures we have considered a simplified example of a centralized 
network with a single center. Such system produces many equilibria, and this dynamical structure can be 
protected against large space dependent, random perturbations. We show that in general, designing an 
optimal network that is protected against such perturbations boils down to finding the minimum energy 
of a spin glass hamiltonian, which is a computationally hard problem. However, for a large number of 
satellites, the randomness is filtered and reliable protection against perturbations results as a solution to 
a quadratic programming problem, that can be solved in polynomial time. We expect that similar results 
hold more generally, for networks with any number of centers. This suggests an evolutionary bias towards 
centralized networks where hubs are subjected to control from many satellites. 

These findings can be interpreted in terms of gene networks. The fiexibility control by satellites, and 
not by transcription factors (centers) can be a major property of such networks. It may be easier to act on 
a satellite (by silencing or reactivating it), then to perform similar actions on a center (deletion of a hub 
proves most of the time to be lethal) . We have proposed miRNAs and CREMs as possible candidates for 
satellite nodes in gene networks controlling pattering in development. A few examples of such centralized 
motifs are known, such for instance the enhancer system of the even-skipped gene of Drosophila (Ludwig 
et al 2011). The process of reconstruction of such networks is only at the beginning (see for instance 
(Berezikov 2011)). One could expect that many more examples of centralized motifs and networks will be 
found during this process. 
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Appendix: Proofs and estimates 

I. The proof of Proposition 2.1 

To outline the proof, let us notice that our system has a typical form, where slow (v) and fast (it) 
components are separated: 

vt — kF{v,u), ut^Au + KG{v). (78) 

Let us present u as u = U + u, where U = —kA^^G{v) and u is a new unknown. Let us notice that Ui 
are solutions of ([T5)) under boundary conditions ([5]) and that \Ui\ < ck. 
By substituting u = U + u into (|78p . we obtain 

vt ^ kF{v,U + u), ut^ Au + k'^Gi{v,U + u), (79) 

where Gi = K^^A^^G'{v)vt = A^^F{v,U + u), the operator A = diag{diA — Xi}. Let us show that 
Gi{u,v) is a uniformly bounded map in H for all u,v satisfying a priori estimates PH)) . For sufficiently 
smooth initial data </), ^ S these estimates and evolution equation © imply 

lk(i)IU <Ci, t>0, ae (0,1). (80) 

The Sobolcv embedding gives then 

l|Vi;(t)|U,(,2)<c||«(t)|U<C2, t>0, ae (1/2,1). (81) 

To estimate now Gi — {wi, ...,wnY^, let us notice that Wi satisfy the following equations: 

{d,A ~ \)w, = g,(a;, v){d^A ~ A,)w,, (82) 

where gi are smooth functions with uniformly bounded derivatives. Our goal is, thus, to estimate ||Vw|| 
through ||Vw||lj and HwHa- Let us multiply (|82p through wi and then integrate the left hand and the 
right hand sides of the obtained equations by parts. We find 

||Vi«,|p < ci||V«;,||||Vw,|| +C2((Vw)M?i;|). (83) 

To estimate ((V^;)^ we use the Cauchy-Schwartz inequality 

|((Vi;)Mu;|)|<c||VH|L,(f2)lkll, 
Now we can apply (j8ip and the Cauchy inequality with a parameter a > that gives 

llVwII ||Vw|| < cia||Vw|p +Ca"i||Vw|p, (84) 
and if a > is small enough {c\a < 1), wc obtain, by ((83|) and ((84| . the need estimate: 

||Vw|| < C3. 

The second equation in ([75]) then entails 

llwllt < -P\\u\\+ K^svcpWGiW, 
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where /3 = min{Ai} > is independent of k. This gives 

\\u(t)\\ < ll"(0)|| cxp(-/3<) + CiK^. 
In a similar way one can obtain the same estimate for ||V?i||. This completes the proof. 

II. Estimates for network viability via spin hamiltonian 

Assmne that for some ^{x) = {£,i{x), ...,^n{x)) there holds 

li{x,a-))<S. (85) 

Let us obtain estimates of deviations v = v ~ q and Ui = Ui — Ui{q), where v = q{x), u = Ui define an 
equilibrium stationary solution for £,i{x) = 0. These estimates hold only due to the special structure of 
our network: we admit that are not small, nonetheless, the summarized effect of these perturbations is 
small. We assume 

Ui(x,0) = 0, u(x,0) = 0. (86) 

Let us present the functions Ui as sums Ui = fji + wt, where fji are defined by (pT|) . For Wi,v we then 
obtain 

= doAd - \ov + a{p{U + u{t)) -h)- a{p{U) - h), (87) 
Biii(f] 

^^^^d,Ai~Kv + FMr),OdT, (88) 
Fi{v, = '^{biiq + v)-h,+ ^i) - u{hq - h, + ^i)- 

11^.11 < c||w||. (89) 



where 

Let us observe that 
Condition ([55]) implies that 



Mri{-))\\<5. (90) 
Let us introduce \\w\\, by WwW^ = X^iLi ll^ilP ^'^^ \'A by |ap = '}2!i=i Wii^ ■ Then 

\\p(w)\\ < \a\\\w\\. 

By dST]), dMl), dMl) and ^ now one obtains inequalities for ||?;||, ||w,H: 

<-Ao||^|P + C3(||p(fy)|| + |a|||^z;||), (91) 



2dt 



d\\ 



w 



|2 



2dt 



< -A||u;|p +C4||w||. (92) 



where min^ = A > 0. Assume that min^ Ai > are large enough. Moreover, for large N the coefficient 
c^\a\ < 1. Combining (^1]) . (|92p one obtains the inequality for = + H^lp: 

< -M\y\? + (Ao - A)||u;|p + C3(5||«|| + csl^H \\v\\. (93) 

We apply now the Cauchy inequality xy < ax^ + a~^y^ to the two terms in the right hand side of this 
last inequality. This gives 

"^II^IP ^ \ in^i|2 I /\ \ML,,I|2 , „-l||„,,||2 I „ „IL--,I|2 , „ „-lr2||--,| 



2dt 



< -Ao||y|r + (Ao - A)||u;|r + a" IkH + ceaM' + C7a-'d^\\v\\. (94) 
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We adjust an a such that cqq < \o/2. If A is large enough, gives then 

\\Yit)\\ <c(5+||r(0)||exp(-Aoi/2). (95) 
Then (pSj) imphes that for large t there holds 

sup < cgS 

t>0 

with a constant cg > 0. This gives us the need estimate of v via the spin hamiltonian. 
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